import pandas as pd
import numpy as np  # 矩阵运算
import matplotlib.pyplot as plt
import seaborn as sns  # 图形显示

plt.style.use('fivethirtyeight')  # 样式美化
from sklearn.metrics import classification_report  # 这个包式评价报告

# 建立一个逻辑回归模型来预测一个学生是否被大学录取，根据两次考试的结果来决定每个人的录取机会
# pyehon实现逻辑回归 目标：建立分类器(求解出三个参数 θ0，θ1，θ2)，即得出分界线，备注 θ1对应exam1 θ2对应exam2
# 设定阈值，根据阈值判断录取结果
# 备注：阈值指的是最终得到的概率值，将概率值转化为一个类别，一般是>0.5被录取，<0.5不录取

# sigmoid 函数：映射到概率的函数
# model：返回预测结果值
# cost：根据参数计算损失
# gradient：计算每个参数的梯度方向
# descent：进行参数更新
# accurancy：计算精度

data = pd.read_csv("ex2data1.txt", names=['exam1', 'exam2', 'admitted'])
print(data.head())
#  数据如下
#        exam1      exam2  admitted
# 0  34.623660  78.024693          0
# 1  30.286711  43.894998          0
# 2  35.847409  72.902198          0
# 3  60.182599  86.308552          1
# 4  79.032736  75.344376          1

print(data.describe())
#             exam1       exam2   admitted
# count  100.000000  100.000000  100.000000
# mean    65.644274   66.221998    0.600000
# std     19.458222   18.582783    0.492366
# min     30.058822   30.603263    0.000000
# 25%     50.919511   48.179205    0.000000
# 50%     67.032988   67.682381    1.000000
# 75%     80.212529   79.360605    1.000000
# max     99.827858   98.869436    1.000000

sns.set(context="notebook", style="darkgrid", palette=sns.color_palette("RdBu", 2))  # 设置样式参数
sns.lmplot(x='exam1', y='exam2', hue='admitted', data=data, height=6, fit_reg=False, scatter_kws={"s": 50})

plt.show()


def get_X(df):  # 读取特征
    ones = pd.DataFrame({"ones": np.ones(len(df))})
    data = pd.concat([ones, df], axis=1)  # 合并数据，根据列合并      axis=1的时候，concat就是行对齐，然后将不同列名称的两张表合并，加列
    # 其实就是相当于在上面的数据前加了一列全1
    return data.iloc[:, :-1].values   # 这个操作返回ndarray，不是矩阵
    # 在GH18458这个补丁中，已经把DataFrame对象和Series对象的as_matrix方法都删除了，使用的时候自然就报错了。
    # 所以使用 .values来代替


def get_Y(df):  # 读取标签
    return np.array(df.iloc[:, -1])  # df.iloc[:, -1] 是指df的最后一列


# 归一化
def normalize_feature(df):
    return df.apply(lambda column: (column - column.mean()) / column.std())  # 特征缩放再逻辑回归同样适用


X = get_X(data)
print(X.shape)
# (100, 3)

y = get_Y(data)
print(y.shape)
# (100,)


# Sigmoid函数
# g 代表一个常用的逻辑函数为S形函数，公式为 g(z) = 1 / (1+e^-z)
# 合起来，我们得到逻辑回归模型的假设函数：hθ(x) = 1 / (1+e^-θTx)


def sigmoid(z):
    return 1 / (1 + np.exp(-z))


fig, ax = plt.subplots(figsize=(8, 6))
ax.plot(np.arange(-10, 10, step=0.01), sigmoid(np.arange(-10, 10, step=0.01)))
ax.set_ylim((-0.1, 1.1))  # lim：轴线显示长度
ax.set_xlabel('z', fontsize=18)
ax.set_ylabel('g(z)', fontsize=18)
ax.set_title('sigmoid function', fontsize=18)
plt.show()


# 代价函数
theta = np.zeros(3)
print(theta)
# [0. 0. 0.]


def cost(theta, X, y):
    return np.mean(-y * np.log(sigmoid(X @ theta)) - (1 - y) * np.log(1 - sigmoid(X @ theta)))


print(cost(theta, X, y))
# 0.6931471805599453


# 梯度下降
def gradient(theta, X, y):
    return (1 / len(X)) * X.T @ (sigmoid(X @ theta) - y)


print(gradient(theta, X, y))
# [ -0.1        -12.00921659 -11.26284221]


# 拟合函数
import scipy.optimize as opt
res = opt.minimize(fun=cost, x0=theta, args=(X, y), method='Newton-CG', jac=gradient)
print(res)
#      fun: 0.20349770451259855
#      jac: array([1.62947970e-05, 1.11339134e-03, 1.07609314e-03])
#  message: 'Optimization terminated successfully.'
#     nfev: 71
#     nhev: 0
#      nit: 28
#     njev: 187
#   status: 0
#  success: True
#        x: array([-25.16576744,   0.20626712,   0.20150754])
#
# Process finished with exit code 0


def predict(x, theta):
    prob = sigmoid(x @ theta)
    return (prob >= 0.5).astype(int)  # 实现变量类型转换


final_theta = res.x
y_pred = predict(X, final_theta)
print(classification_report(y, y_pred))
#               precision    recall  f1-score   support
#
#            0       0.87      0.85      0.86        40
#            1       0.90      0.92      0.91        60
#
#     accuracy                           0.89       100
#    macro avg       0.89      0.88      0.88       100
# weighted avg       0.89      0.89      0.89       100


# 寻找决策边界
print(res.x)
# [-25.16576744   0.20626712   0.20150754]

coef = -(res.x / res.x[2])
print(coef)
# [124.88747248  -1.02361985  -1.        ]

x = np.arange(130, step=0.1)
y = coef[0] + coef[1]*x


# 绘制决策边界
sns.set(context="notebook", style="ticks", font_scale=1.5)

sns.lmplot(x='exam1', y='exam2', hue='admitted', data=data, height=6, fit_reg=False, scatter_kws={"s": 25})
plt.plot(x, y, 'grey')
plt.xlim(0, 130)
plt.ylim(0, 130)
plt.title('Decision Boundary')
plt.show()

